setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
gwas_fn=list.files(pattern = "txt",path = "E:/0 公共数据库差异情况/db_for_5hmc/MAGMA分析的结果",full.names=T)
con.file=list.files(pattern = "hg19_multianno.csv",path = "E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207",full.names=T)

file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

caseid=unique(unlist(strsplit(file3$Gene.refGene,";")))
caseid=caseid[!caseid=="NONE"]

n=seq(-6,-2,1)
threshold=as.numeric(c(10^n,5*10^n))###设置几个阈值

col_names=c("threshold","GWAS_MAGMA","con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=9))
names(result)=col_names
result=result[-1,]

for(l in 1:length(threshold)){
for(j in 1:length(gwas_fn)){
gwas.file=read.table(gwas_fn[j],head=T,sep="\t")
gwas.file_sig=gwas.file[gwas.file$P<threshold[l],]
sig_id=unique(gwas.file_sig$V6)

for(k in 1:length(con.file)){
  con=read.csv(con.file[k],header=T)
  conid=unique(unlist(strsplit(con$Gene.refGene,";")))
  conid=conid[!conid=="NONE"]
  result_tmp=data.frame(matrix(NA,1,ncol=9))
  names(result_tmp)=col_names
  result_tmp[1,1]=threshold[l]
  result_tmp[1,2]=gsub(".result.txt","",gsub("E:/0 公共数据库差异情况/db_for_5hmc/MAGMA分析的结果/","",gwas_fn[j]))
  result_tmp[1,3]=gsub("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207/","",gsub(".anno.hg19_multianno.csv","",con.file[k]))
  result_tmp[1,4]=length(intersect(caseid,sig_id))
  result_tmp[1,5]=length(caseid)-length(intersect(caseid,sig_id))
  result_tmp[1,6]=length(intersect(conid,sig_id))
  result_tmp[1,7]=length(conid)-length(intersect(conid,sig_id))
  result_tmp[1,8]=fisher.test(matrix(c(result_tmp[1,4],result_tmp[1,5],result_tmp[1,6],result_tmp[1,7]),nrow = 2))$estimate
  result_tmp[1,9]=fisher.test(matrix(c(result_tmp[1,4],result_tmp[1,5],result_tmp[1,6],result_tmp[1,7]),nrow = 2))$p.value
  result=rbind(result,result_tmp)
}
}
}
write.csv(result,"20201209/GWAS_MAGMA_enrichment.csv",quote = F,row.names = F)

library(meta)
col_names=c("threshold","GWAS_MAGMA","case_overlap","OR","p.value","upper","lower")
rt=data.frame(matrix(NA,1,ncol=7))
names(rt)=col_names
rt=rt[-1,]

for(i in seq(1,90,3)){
data.tmp=data.frame(result[i:(i+2),])
data.tmp$case_not=data.tmp$case_overlap+data.tmp$case_not
data.tmp$con_not=data.tmp$con_overlap+data.tmp$con_not
metaor3<-metabin(case_overlap,case_not,con_overlap,con_not,data=data.tmp,sm="OR")

rt.tmp=data.frame(matrix(NA,1,ncol=7))
names(rt.tmp)=col_names
rt.tmp[1,1]=data.tmp[1,1]
rt.tmp[1,2]=data.tmp[1,2]
rt.tmp[1,3]=data.tmp[1,4]
rt.tmp[1,4]=OR=exp(metaor3$TE.fixed)
rt.tmp[1,5]=metaor3$pval.fixed
rt.tmp[1,6]=upper=exp(metaor3$upper.fixed)
rt.tmp[1,7]=lower=exp(metaor3$lower.fixed)
rt=rbind(rt,rt.tmp)
}
rt$FDR=p.adjust(rt$p.value,method = "BH")
write.csv(rt,"20201209/meta.p.GWAS_MAGMA_enrichment.csv",quote = F,row.names = F)